home *** CD-ROM | disk | FTP | other *** search
/ Computer Select (Limited Edition) / Computer Select.iso / dobbs / v17n04 / wavelet.exe / WAVEDEMO.C < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-17  |  15.5 KB  |  626 lines
  1. #define WAVEDEMO
  2.  
  3. #include <conio.h>
  4. #include <math.h>
  5. #include <stdlib.h>
  6. #include <stdio.h>
  7. #include <bios.h>
  8. #include <dos.h>
  9.  
  10. typedef enum { FALSE, TRUE } tftype;
  11. typedef enum { APPROX, DETAIL, INPUT, OUTPUT } signaltype;
  12. typedef enum { DECOMP, RECON } wavetype;
  13. typedef enum { NONE, HELP1, HELP2 } helptype;
  14.  
  15. #include "wavelet.h"
  16. #include "wave_mgt.h"
  17. #include "wavegrph.h"
  18.  
  19. void SignalFileIn(char *filename, int *length, int levels);
  20. void SignalFileOut(char *filename, int length);
  21. void CoeffFileIn(char *filename, int *length, int *levels,
  22.                                                                                     double *alpha, double *beta);
  23. void CoeffFileOut(char *filename, int signal_length, int levels,
  24.                                                                                     double alpha, double beta);
  25.  
  26. double InData[317], OutData[317];
  27. double alphaD, alphaR, betaD, betaR;
  28. double g[6], gtilda[6], h[6], htilda[6];
  29. double WaveCoeffD[6], WaveCoeffR[6];
  30. double **ExCoeffD, **ExCoeffR;
  31. int WaveLevels, DataLength, filtlenD, filtlenR;
  32. char InSignalname[20], OutSignalname[20];
  33. char InCoeffname[20], OutCoeffname[20];
  34.  
  35. void SignalFileIn(char *filename, int *length, int levels)
  36. {
  37.     int i, new_length, str_len;
  38.     char buffer[15];
  39.     FILE *in;
  40.  
  41.     str_len = GetGraphicInput(384, trow1, buffer, 14);
  42.     if (str_len && ((in = fopen(buffer, "rt")) == NULL))
  43.     {
  44.         GraphMessage(328, trow7, "CANNOT OPEN THIS FILE         ");
  45.         beep();
  46.         delay(1000);
  47.     }
  48.  
  49.     if (in && str_len)
  50.     {
  51.         i = 0;
  52.         do  /* assign the new file to the input data */
  53.         {
  54.             filename[i] = buffer[i];
  55.         } while (buffer[i++] != '\0');
  56.  
  57.         fscanf(in, "%d", &new_length); /* get number of data points */
  58.         if (new_length != *length)  /* resize dynamic storage if necessary */
  59.         {
  60.             DestroyTreeStorage(ExCoeffR, levels);
  61.             DestroyTreeStorage(ExCoeffD, levels);
  62.             ExCoeffD = BuildTreeStorage(new_length, levels);
  63.             ExCoeffR = BuildTreeStorage(new_length, levels);
  64.             *length = new_length;
  65.         }
  66.  
  67.         for (i = 0; i < *length; i++) /* load in the data points */
  68.             fscanf(in, "%lf", &InData[i]);
  69.  
  70.         fclose(in);
  71.         for (i = *length; i < 312; i++)
  72.             InData[i] = 0.0;
  73.  
  74.     }
  75. }
  76.  
  77.  
  78. void SignalFileOut(char *filename, int length)
  79. {
  80.     int i, str_len;
  81.     char buffer[15];
  82.     FILE *out;
  83.  
  84.     str_len = GetGraphicInput(392, trow1, buffer, 14);
  85.     if (str_len && ((out = fopen(buffer, "wt")) == NULL))
  86.     {
  87.         GraphMessage(328, trow7, "CANNOT OPEN THIS FILE         ");
  88.         beep();
  89.         delay(1000);
  90.     }
  91.  
  92.     if (out && str_len)
  93.     {
  94.         i = 0;
  95.         do /* assign the new file to the output data */
  96.         {
  97.             filename[i] = buffer[i];
  98.         } while (buffer[i++] != '\0');
  99.  
  100.         fprintf(out, "%d\n", length); /* output number of data points */
  101.         for (i = 0; i < length; i++)
  102.             fprintf(out, "%20.14le\n", OutData[i]); /* output the data */
  103.  
  104.         fclose(out);
  105.     }
  106. }
  107.  
  108.  
  109. void CoeffFileIn(char *filename, int *length, int *levels,
  110.                                                                                     double *alpha, double *beta)
  111. {
  112.     int i, j, new_length, new_levels, str_len;
  113.     char buffer[15];
  114.     FILE *in;
  115.  
  116.     str_len = GetGraphicInput(384, trow3, buffer, 14);
  117.     if (str_len && ((in = fopen(buffer, "rt")) == NULL))
  118.     {
  119.         GraphMessage(328, trow7, "CANNOT OPEN THIS FILE         ");
  120.         beep();
  121.         delay(1000);
  122.     }
  123.  
  124.     if (in && str_len)
  125.     {
  126.         i = 0;
  127.         do /* assign the new file to the input coefficients */
  128.         {
  129.             filename[i] = buffer[i];
  130.         } while (buffer[i++] != '\0');
  131.  
  132.          /* get length of transform and number of levels */
  133.         fscanf(in, "%d   %d\n", &new_length, &new_levels);
  134.          /* resize dynamic storage if necessary */
  135.         if ((new_length != *length) || (new_levels !=  *levels))
  136.         {
  137.             DestroyTreeStorage(ExCoeffR, *levels);
  138.             DestroyTreeStorage(ExCoeffD, *levels);
  139.             ExCoeffD = BuildTreeStorage(new_length, new_levels);
  140.             ExCoeffR = BuildTreeStorage(new_length, new_levels);
  141.             *length = new_length;
  142.             *levels = new_levels;
  143.         }
  144.  
  145.         fscanf(in, "%le %le", alpha, beta); /* load values of alpha and beta */
  146.  
  147.         for (i = 0; i < *levels; i++) /* load detail coefficients */
  148.         {
  149.             new_length = (*length >> (i + 1)) + 5;
  150.             for (j = 0; j < new_length; j++)
  151.                 fscanf(in, "%le", &ExCoeffR[(2 * i) + 1][j]);
  152.         }
  153.  
  154.         for (j = 0; j < new_length; j++) /* load approximation coefficients */
  155.             fscanf(in, "%le", &ExCoeffR[2 * (*levels - 1)][j]);
  156.  
  157.         fclose(in);
  158.     }
  159. }
  160.  
  161.  
  162. void CoeffFileOut(char *filename, int signal_length, int levels,
  163.                                                                                     double alpha, double beta)
  164. {
  165.     int i, j, length, str_len;
  166.     char buffer[15];
  167.     FILE *out;
  168.  
  169.     str_len = GetGraphicInput(392, trow3, buffer, 14);
  170.     if (str_len && ((out = fopen(buffer, "wt")) == NULL))
  171.     {
  172.         GraphMessage(328, trow7, "CANNOT OPEN THIS FILE         ");
  173.         beep();
  174.         delay(1000);
  175.     }
  176.  
  177.     if (out && str_len)
  178.     {
  179.         i = 0;
  180.         do /* assign the new file to the output coefficients */
  181.         {
  182.             filename[i] = buffer[i];
  183.         } while (buffer[i++] != '\0');
  184.  
  185.             /* write length of transform and number of levels */
  186.         fprintf(out, "%d   %d\n", signal_length, levels);
  187.             /* write values of alpha and beta */
  188.         fprintf(out, "%20.14le   %20.14le\n", alpha, beta);
  189.  
  190.         for (i = 0; i < levels; i++) /* write detail coefficients */
  191.         {
  192.             length = (signal_length >> (i + 1)) + 5;
  193.             for (j = 0; j < length; j++)
  194.                 fprintf(out, "%20.14le\n", ExCoeffD[(2 * i) + 1][j]);
  195.         }
  196.  
  197.         for (j = 0; j < length; j++) /* write approximation coefficients */
  198.             fprintf(out, "%20.14le\n", ExCoeffD[2 * (levels - 1)][j]);
  199.  
  200.         fclose(out);
  201.     }
  202. }
  203.  
  204.  
  205. void main(void)
  206. {
  207.     tftype quit;
  208.     double temp, pi = 3.14159265358979323846;
  209.     int i;
  210.     char buffer[20];
  211.     unsigned char str_len;
  212.  
  213.     union
  214.     {
  215.         char c[2];
  216.         int i;
  217.     } key;
  218.  
  219.     if (GraphSetup())
  220.         exit(255);
  221.  
  222.     DataLength = 312;
  223.     DispLength = DataLength;
  224.     for (i = 0; i < DataLength; i++)
  225.         InData[i] = sin(0.02 * pi * (double) i);
  226.  
  227.     TopResShown = 0;
  228.     WaveLevels = 3;
  229.     DispWaveLevels = WaveLevels;
  230.     ExCoeffD = BuildTreeStorage(DataLength, WaveLevels);
  231.     ExCoeffR = BuildTreeStorage(DataLength, WaveLevels);
  232.     DispExCoeff = ExCoeffD;
  233.     DispCoeffname = OutCoeffname;
  234.     DispData = InData;
  235.     DispSignalname = InSignalname;
  236.     DispAlpha = &alphaD;
  237.     DispBeta = &betaD;
  238.     DispWaveCoeff = WaveCoeffD;
  239.     WaveShown = DECOMP;
  240.     HelpToggle = HELP1;
  241.     WaveletCoeffs(alphaD, betaD, WaveCoeffD);
  242.     filtlenD = MakeWaveletFilters(WaveCoeffD, htilda, gtilda, DECOMP);
  243.     WaveletCoeffs(alphaR, betaR, WaveCoeffR);
  244.     filtlenR = MakeWaveletFilters(WaveCoeffR, h, g, RECON);
  245.     DrawViewedLevels();
  246.  
  247.     quit = FALSE;
  248.     delay(0);
  249.  
  250.     do
  251.     {
  252.         if (bioskey(1))
  253.         {
  254.             key.i = bioskey(0);
  255.             if (!key.c[0])
  256.             {
  257.                 switch (key.c[1])
  258.                 {
  259.                     case 0x48:    /* Up Arrow */
  260.                         if (ViewingData && (ViewStart > 0))
  261.                             ShowViewData(-1);
  262.                         else
  263.                             beep();
  264.  
  265.                         break;
  266.                     case 0x49:    /* PGUP */
  267.                         if (!ViewingData && TopResShown)
  268.                             ScrollInLevel(--TopResShown, WaveLevels, 1);
  269.                         else if ((ViewingData == TRUE) && (ViewStart > 0))
  270.                             ShowViewData(-(viewlines - 1));
  271.                         else
  272.                             beep();
  273.  
  274.                         break;
  275.                     case 0x50:    /* Dn Arrow */
  276.                         if (ViewingData && (ViewStart < (ViewLength - 1)))
  277.                             ShowViewData(1);
  278.                         else
  279.                             beep();
  280.  
  281.                         break;
  282.                     case 0x51:    /* PGDN */
  283.                         if (!ViewingData && (TopResShown < WaveLevels))
  284.                             ScrollInLevel(3 + TopResShown++, WaveLevels, 0);
  285.                         else if ((ViewingData == TRUE)
  286.                             && (ViewStart < (ViewLength - 1)))
  287.                             ShowViewData(viewlines - 1);
  288.                         else
  289.                             beep();
  290.  
  291.                         break;
  292.                     case 0x3B:    /* F1 - HELP */
  293.                         if ((TopResShown == 0) && (HelpToggle != NONE) && !ViewingData)
  294.                         {
  295.                             HelpToggle = NONE;
  296.                             DrawPanel();
  297.                         }
  298.                         else if ((TopResShown == 0) && !HelpToggle && !ViewingData)
  299.                         {
  300.                             HelpToggle = HELP1;
  301.                             DrawPanel();
  302.                         }
  303.                         else
  304.                             beep();
  305.  
  306.                         break;
  307.                     default:
  308.                         beep();
  309.                         delay(50);
  310.                         beep();
  311.                 }
  312.             }
  313.             else
  314.             {
  315.                 switch (key.c[0])
  316.                 {
  317.                     case 0x1b:    /* ESC exit the program */
  318.                         if (ViewingData)
  319.                         {
  320.                             ExitViewData();
  321.                             DrawViewedLevels();
  322.                         }
  323.                         else
  324.                             quit = TRUE;
  325.  
  326.                         break;
  327.                     case 'A': /* set value of alpha */
  328.                     case 'a':
  329.                         if ((TopResShown == 0) && !HelpToggle && !ViewingData)
  330.                         {
  331.                             str_len = GetGraphicInput(376, trow4, buffer, 15);
  332.                             temp = atof(buffer);
  333.                             if ((WaveShown == DECOMP) && str_len)
  334.                             {
  335.                                 if ((temp >= -pi) && (temp <= pi))
  336.                                 {
  337.                                     alphaD = temp;
  338.                                     WaveletCoeffs(alphaD, betaD, WaveCoeffD);
  339.                                     filtlenD = MakeWaveletFilters(WaveCoeffD,
  340.                                                                             htilda, gtilda, DECOMP);
  341.                                 }
  342.                                 else
  343.                                     beep();
  344.                             }
  345.                             else if (str_len)
  346.                             {
  347.                                 if ((temp >= -pi) && (temp <= pi))
  348.                                 {
  349.                                     alphaR = temp;
  350.                                     WaveletCoeffs(alphaR, betaR, WaveCoeffR);
  351.                                     filtlenR = MakeWaveletFilters(WaveCoeffR, h, g, RECON);
  352.                                 }
  353.                                 else
  354.                                     beep();
  355.                             }
  356.  
  357.                             DrawPanel();
  358.                         }
  359.                         else
  360.                             beep();
  361.  
  362.                         break;
  363.                     case 'B': /* set value of beta */
  364.                     case 'b':
  365.                         if ((TopResShown == 0) && !HelpToggle && !ViewingData)
  366.                         {
  367.                             str_len = GetGraphicInput(368, trow5, buffer, 15);
  368.                             temp = atof(buffer);
  369.                             if ((WaveShown == DECOMP) && str_len)
  370.                             {
  371.                                 if ((temp >= -pi) && (temp <= pi))
  372.                                 {
  373.                                     betaD = temp;
  374.                                     WaveletCoeffs(alphaD, betaD, WaveCoeffD);
  375.                                     filtlenD = MakeWaveletFilters(WaveCoeffD,
  376.                                                                             htilda, gtilda, DECOMP);
  377.                                 }
  378.                                 else
  379.                                     beep();
  380.                             }
  381.                             else if (str_len)
  382.                             {
  383.                                 if ((temp >= -pi) && (temp <= pi))
  384.                                 {
  385.                                     betaR = temp;
  386.                                     WaveletCoeffs(alphaR, betaR, WaveCoeffR);
  387.                                     filtlenR = MakeWaveletFilters(WaveCoeffR, h, g, RECON);
  388.                                 }
  389.                                 else
  390.                                     beep();
  391.                             }
  392.  
  393.                             DrawPanel();
  394.                         }
  395.                         else
  396.                             beep();
  397.  
  398.                         break;
  399.                     case 'C': /* copy expansion coefficients */
  400.                     case 'c':
  401.                         if (!HelpToggle && (WaveShown == DECOMP) && !ViewingData)
  402.                         {
  403.                             TreeCopy(ExCoeffD, ExCoeffR, DataLength, WaveLevels);
  404.                             DrawViewedLevels();
  405.                         }
  406.                         else if (!HelpToggle
  407.                             && (WaveShown == RECON) && !ViewingData)
  408.                         {
  409.                             TreeCopy(ExCoeffR, ExCoeffD, DataLength, WaveLevels);
  410.                             DrawViewedLevels();
  411.                         }
  412.                         else
  413.                             beep();
  414.  
  415.                         break;
  416.                     case 'E': /* execute displayed transform */
  417.                     case 'e':
  418.                         if (!HelpToggle && (WaveShown == DECOMP) && !ViewingData)
  419.                         {
  420.                             WaveletDecomposition(InData, DataLength, htilda,
  421.                                                     gtilda, filtlenD, WaveLevels, ExCoeffD);
  422.                             DispMSE = CalculateMSE(InData, OutData, DataLength);
  423.                             DrawViewedLevels();
  424.                         }
  425.                         else if (!HelpToggle
  426.                 && (WaveShown == RECON) && !ViewingData)
  427.                         {
  428.                             WaveletReconstruction(ExCoeffR, DataLength, h,
  429.                                                         g, filtlenR, WaveLevels, OutData);
  430.                             DispMSE = CalculateMSE(InData, OutData, DataLength);
  431.                             DrawViewedLevels();
  432.                         }
  433.                         else
  434.                             beep();
  435.  
  436.                         break;
  437.                     case 'I': /* input signal read */
  438.                     case 'i':
  439.                         if (!HelpToggle && (TopResShown == 0)
  440.                                 && (WaveShown == DECOMP) && !ViewingData)
  441.                         {
  442.                             SignalFileIn(InSignalname, &DataLength, WaveLevels);
  443.                             DispLength = DataLength;
  444.                             DispExCoeff = ExCoeffD;
  445.                             DrawViewedLevels();
  446.                         }
  447.                     else
  448.                             beep();
  449.  
  450.                         break;
  451.                     case 'L': /* set number of levels in transform */
  452.                     case 'l':
  453.                         if ((TopResShown == 0) && !HelpToggle && !ViewingData)
  454.                         {
  455.                             str_len = GetGraphicInput(384, trow6, buffer, 1);
  456.                             i = atoi(buffer);
  457.                             if ((i > 0) && (i < 7) && str_len)
  458.                             {
  459.                                 DestroyTreeStorage(ExCoeffR, WaveLevels);
  460.                                 DestroyTreeStorage(ExCoeffD, WaveLevels);
  461.                                 WaveLevels = i;
  462.                                 DispWaveLevels = WaveLevels;
  463.                                 ExCoeffD = BuildTreeStorage(DataLength, WaveLevels);
  464.                                 ExCoeffR = BuildTreeStorage(DataLength, WaveLevels);
  465.                                 if (WaveShown == DECOMP)
  466.                                     DispExCoeff = ExCoeffD;
  467.                                 else
  468.                                     DispExCoeff = ExCoeffR;
  469.                             }
  470.                             else
  471.                                 beep();
  472.  
  473.                             DrawViewedLevels();
  474.                         }
  475.                         else
  476.                             beep();
  477.  
  478.                         break;
  479.                     case 'M': /* show more help */
  480.                     case 'm':
  481.                         if ((TopResShown == 0) && (HelpToggle == HELP1) && !ViewingData)
  482.                         {
  483.                             HelpToggle = HELP2;
  484.                             DrawPanel();
  485.                         }
  486.                         else if ((TopResShown == 0) && (HelpToggle == HELP2))
  487.                         {
  488.                             HelpToggle = HELP1;
  489.                             DrawPanel();
  490.                         }
  491.                         else
  492.                             beep();
  493.  
  494.                         break;
  495.                     case 'O': /* output signal write */
  496.                     case 'o':
  497.                         if (!HelpToggle && (TopResShown == 0)
  498.                 && (WaveShown == RECON) && !ViewingData)
  499.                         {
  500.                             SignalFileOut(OutSignalname, DataLength);
  501.                             DrawPanel();
  502.                         }
  503.                         else
  504.                             beep();
  505.  
  506.                         break;
  507.                     case 'P': /* prune selected levels of expansion coeffs from tree */
  508.                     case 'p':
  509.                         if (!HelpToggle && (TopResShown == 0)
  510.                 && (WaveShown == RECON) && !ViewingData)
  511.                         {
  512.                             GraphMessage(424, trow6, "ZERO DETAIL IN TOP   LEVELS");
  513.                             str_len = GetGraphicInput(550, trow6, buffer, 1);
  514.                             i = atoi(buffer);
  515.                             if ((i > 0) && (i <= WaveLevels) && str_len)
  516.                                 ZeroTreeDetail(DispExCoeff, DataLength, i);
  517.  
  518.                             DrawViewedLevels();
  519.                         }
  520.                         else
  521.                             beep();
  522.  
  523.                         break;
  524.                     case 'R': /* read expansion coefficients */
  525.                     case 'r':
  526.                         if ((TopResShown == 0) && !HelpToggle
  527.                 && (WaveShown == RECON) && !ViewingData)
  528.             {
  529.                             CoeffFileIn(InCoeffname, &DataLength, &WaveLevels,
  530.                                                                                                     &alphaR, &betaR);
  531.                             WaveletCoeffs(alphaR, betaR, WaveCoeffR);
  532.                             filtlenR = MakeWaveletFilters(WaveCoeffR, h, g, RECON);
  533.                             DispLength = DataLength;
  534.                             DispWaveLevels = WaveLevels;
  535.                             DispExCoeff = ExCoeffR;
  536.                             DrawViewedLevels();
  537.                         }
  538.                         else
  539.                             beep();
  540.  
  541.                         break;
  542.                     case 'T':
  543.                     case 't':
  544.                         if (!ViewingData)
  545.                         {
  546.                             if (WaveShown == DECOMP)
  547.                             {
  548.                                 DispData = OutData;
  549.                                 DispSignalname = OutSignalname;
  550.                                 DispAlpha = &alphaR;
  551.                                 DispBeta = &betaR;
  552.                                 DispWaveCoeff = WaveCoeffR;
  553.                                 DispExCoeff = ExCoeffR;
  554.                                 DispCoeffname = InCoeffname;
  555.                                 WaveShown = RECON;
  556.                             }
  557.                             else
  558.                             {
  559.                                 DispData = InData;
  560.                                 DispSignalname = InSignalname;
  561.                                 DispAlpha = &alphaD;
  562.                                 DispBeta = &betaD;
  563.                                 DispWaveCoeff = WaveCoeffD;
  564.                                 DispExCoeff = ExCoeffD;
  565.                                 DispCoeffname = OutCoeffname;
  566.                                 WaveShown = DECOMP;
  567.                             }
  568.  
  569.                             DrawViewedLevels();
  570.                         }
  571.                         else
  572.                             beep();
  573.  
  574.                         break;
  575.                     case 'V': /* view the expansion coefficient */
  576.                     case 'v':
  577.                         if (!HelpToggle && !ViewingData)
  578.                         {
  579.                             InitViewData(InData, OutData, ExCoeffD, ExCoeffR);
  580.                             ShowViewData(0);
  581.                         }
  582.                         else
  583.                             beep();
  584.  
  585.                         break;
  586.                     case 'W': /* write expansion coefficients */
  587.                     case 'w':
  588.                         if (!HelpToggle && (TopResShown == 0)
  589.                 && (WaveShown == DECOMP) && !ViewingData)
  590.                         {
  591.                             CoeffFileOut(OutCoeffname, DataLength, WaveLevels,
  592.                                                                                                             alphaD, betaD);
  593.                             DrawPanel();
  594.                         }
  595.                         else
  596.                             beep();
  597.  
  598.                         break;
  599.                     case 'Z': /* zero expansion coefficients shown */
  600.                     case 'z':
  601.                         if (!HelpToggle && (WaveShown == DECOMP) && !ViewingData)
  602.                         {
  603.                             TreeZero(ExCoeffD, DataLength, WaveLevels);
  604.                             DrawViewedLevels();
  605.                         }
  606.                         else if (!HelpToggle && (WaveShown == RECON))
  607.                         {
  608.                             TreeZero(ExCoeffR, DataLength, WaveLevels);
  609.                             DrawViewedLevels();
  610.                         }
  611.                         else
  612.                             beep();
  613.  
  614.                         break;
  615.                     default:
  616.                         beep();
  617.                         delay(50);
  618.                         beep();
  619.                 }
  620.             }
  621.         }
  622.     } while (!quit);
  623.  
  624.     GraphExit();
  625. }
  626.